#include "random_cascade_classifier.h"

#include <limits>
#include <fstream>
#include <math.h>
#include <iostream>

#define EPSILON 1.0e-010

class Layer {
public:
	int l;
	svm_node * w;
	double threshold;
	~Layer();
};

Layer::~Layer() {
	delete [] w;
}

void RandomCascadeClassifier::train(svm_problem & prob, int dim, int layer_num, double g) {
	gamma = g;
	vector<svm_node *> positive_samples;
	get_positive_samples(prob, positive_samples);

	layers.resize(layer_num);
	for(int i = 0; i < layer_num; i++) {
		Layer * layer = new Layer();
		layer->w = get_direction(dim, layer->l);
		layer->threshold = get_threshold(positive_samples, layer->w);
		layers[i] = layer;

		if(i % 200 == 0) {
			cout << i << " " << endl;
		}
	}
}

double RandomCascadeClassifier::test(svm_node * x) {
	int passed = 0;
	int filtered = 0;
	for(unsigned i = 0; i < layers.size(); i++) {
		Layer * layer = layers[i];
		double y = test(x, layer);
		
		int prediction = y >= -EPSILON ? 1 : 0;
		if(prediction == 1) {
			passed ++;
		} else {
			filtered ++;
		}
	}
	double ratio = passed / (double) (passed + filtered);
	return ratio;
}

double RandomCascadeClassifier::test(svm_node * x, Layer * layer) {
	double value = chi_square_kernel(x, layer->w);
	double y = value - layer->threshold;
	return y;
}

void RandomCascadeClassifier::dump(char * filename) {
	FILE *output = fopen(filename,"w");

	fprintf(output,"%d %.17g\n",layers.size(), gamma);
	for(unsigned i = 0; i < layers.size(); i++) {
		Layer * layer = layers[i];
		fprintf(output,"%.17g %d", layer->threshold, layer->l);
		svm_node * w = layer->w;
		while(w->index != -1) {
			fprintf(output," %d:%.17g", w->index, w->value);
			w++;
		}
		fprintf(output,"\n");
	}

	fclose(output);
}


void RandomCascadeClassifier::load(char * filename) {
	ifstream fin(filename);
	if(fin.fail()) {
		cout << "Error! " << filename << " not found!" << endl;
		return;
	}

	int n_layer;
	fin >> n_layer >> gamma;

	layers.resize(n_layer);
	for(int i = 0; i < n_layer; i++) {
		Layer * layer = new Layer();
		
		fin >> layer->threshold >> layer->l;
		layer->w = new svm_node [layer->l + 1];
		for(int j = 0; j < layer->l; j++) {
			char colon;
			fin >> layer->w[j].index >> colon >> layer->w[j].value;
		}
		layer->w[layer->l].index = -1;
		layer->w[layer->l].value = 0;
		layers[i] = layer;
	}

	fin.close();
}

void RandomCascadeClassifier::get_positive_samples(svm_problem & prob, vector<svm_node *> & positive) {
	for(int i = 0; i < prob.l; i++) {
		if(prob.y[i] == 1) {
			positive.push_back(prob.x[i]);
		}
	}
}

svm_node * RandomCascadeClassifier::get_direction(int dim, int & l) {
	vector<int> key_list;
	vector<double> value_list;
	for(int i = 0; i < dim; i++) {
		int raw = rand();
		if(raw != 0) {
			key_list.push_back(i + 1);
			value_list.push_back((raw/(double)RAND_MAX - 0.5) * 2);
		}
	}

	l = (int)key_list.size();

	// normalization
	double sum = 0;
	for(int i = 0; i < l; i++) {
		double value = value_list[i];
		sum += value * value;
	}
	sum = pow(sum, 0.5);

	for(int i = 0; i < l; i++) {
		value_list[i] = value_list[i] / sum;
	}

	// construction
	svm_node *x = new svm_node [l + 1];
	for(int i = 0; i < l; i++) {
		x[i].index = key_list[i];
		x[i].value = value_list[i];
	}
	x[l].index = -1;
	x[l].value = 0;

	return x;
}

double RandomCascadeClassifier::get_threshold(vector<svm_node *> & X, const svm_node * w) {
	double threshold = numeric_limits<double>::max();

	for(unsigned i = 0; i < X.size(); i++) {
		svm_node * x = X[i];
		double value = chi_square_kernel(x, w);
		if(value < threshold) {
			threshold = value;
		}
	}
	return threshold;
}

double RandomCascadeClassifier::chi_square_kernel(const svm_node *px, const svm_node *py) {
	double sum = 0;

	while(px->index != -1 && py->index != -1)
	{
		if(px->index == py->index)
		{
			sum += (px->value - py->value)*(px->value - py->value)/(px->value + py->value);
			++px;
			++py;
		}
		else
		{
			if(px->index > py->index)
			{
				sum += py->value*py->value/py->value;
				++py;
			}
			else
			{
				sum += px->value*px->value/px->value;
				++px;
			}
		}
	}

	while(px->index != -1)
	{
		sum += px->value*px->value/px->value;
		++px;
	}

	while(py->index != -1)
	{
		sum += py->value*py->value/py->value;
		++py;
	}

	sum = sum/2.0;

	double result = exp(-gamma*sum);

	return result;	
}

RandomCascadeClassifier::~RandomCascadeClassifier() {
	for(unsigned i = 0; i < layers.size(); i++) {
		delete layers[i];
	}
}
